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Abstract 

We analysed Chandra observations of the bright Fermi pulsar J0633+0632 and found evidence of an absorp¬ 
tion feature in its spectrum at 804^26 (i^he errors here and below are at 90% confidence) with equivalent 
width of bSlgg eV. In addition, we analysed in detail the X-ray spectral continuum taking into account 
co rrelations betwee n the interstellar absorption and the distance to the source. We confirm early findings 
by I Ray et al.l (|201lh that the spectrum contains non-thermal and thermal components. The latter is equally 
well described by the blackbody and magnetised atmosphere models and can be attributed to the emission 
from the bulk of the stellar surface in both cases. The distance to the pulsar is constrained in a range of 
1-4 kpc from the spectral fits. We infer the blackbody surface temperature of lOSlJ^ eV, while for the at¬ 
mosphere model, the temperature, as seen by a distant observer, is eV. In the latter case J0633+0632 

is one of the coldest middle-aged isolated neutron stars with measured temperatures. Finally, it powers an 
extended pulsar wind nebula whose shape suggests a high pulsar proper motion. Looking backwards the 
direction of the presumed proper motion we found a likely birthplace of the pulsar - the Rosette nebula, 
a 50-Myr-old active star-forming region located at about 1?5 from the pulsar. If true, this constrains the 
distance to the pulsar in the range of 1.2-1.8 kpc. 

Keywords: stars: neutron - pulsars: general - pulsars: individual: PSR J0633+0632. 


1 INTRODUCTION 

Usually, X-ray spectra of isolated neutron stars (NSs) 
are well described by a featureless continuum which con¬ 
tains nonthermal and/or thermal component. In rare 
cases, however, absorption features are detected. Un¬ 
derstanding origins of these features is thought to be 
important for various aspects of the NS physics. For in¬ 
stance, they can result from atomic transitions in th e 
mid-Z element NS atmospheres (e.g., iMori fc Hollinnl . 
In this case, as in ordinary stars, it is possible to mea¬ 
sure the surface gravitational redshift and, hence, the 
stellar mass to radius ratio. This is important for in¬ 
dependent diagnostic of the equation of state (EOS) of 
dense matter inside NSs. Absorption features can be 
also identified with either proton or electron cyclotron 
lines. Elec tron cyclotron lines in NS spectra were pre¬ 
dicted by Gnedin fc Snnvaev ( 19741 ) and then discov¬ 
ered by iTrnemper et ^ ( 19781 ) in the Her X-1 binary 
system. These lines were dete cted in many accreting X- 
ray p ulsars since then, (e.g., iRevnivtsev fc Mereghettil 
20141 ) allowing for direct measurements of NS magnetic 


fields. 

Eor isolated NSs, until recently, absorption fea¬ 
tures have been seen in X-ray spectra of only a few 


atypical, radio-silent, pure thermally emitting sources. 
This includes two low-magnetic-field central compact 
objects (C COs) in supernova remnants (SNRs), IE 
1207-5209 (ISa,nwa,l et, 3,1.112002^ and PSR J0821-4300 
( Gotthelf fc Halnern l2009ll . five objects with larger 
fields from the “Magnificent Seven” family (e.g., 
Fires et al. 2Ql4 and references t herein) and one sof t 


gamma repeater SGR 0418+5729 ( Tiengo et al.l 12013 ). 
Sole exceptio n is an ordinary middle -aged radio pulsar 
J1740+1000 ( Kargaltsev et al.|[2012[) . 

Since the launch of the Fe rmi 7 -ray observ atory, 
several dozens of new pulsars ( Abdo et al I liUl have 
been discovered. A substantial number of them are not 
seen in radio but are identified in X-rays. The radio¬ 
quiet PSR J0633+0632 (hereafter J0633) was discov¬ 
ered in a blind search for p ulsations in the Fermi- 
LAT data ( Abdo et al.l 1200^ . Among Ferm+pulsars, 


J0633 is one of the bright est in X-rays, w ith a flux 
Fx ^ 10“^^ erg cm“^ s“^ ( Rav et ^l201l[) . A pulsar 
period P = 297.4 ms and a rotation frequency deriva¬ 
tive imply a characteristic age r = 59.2 kyr, a spin- 
down luminosity F = 1.2 x 10^^ er g s~^ and a surfac e 
magnetic field B = 4.9 x 10^^ G ([Abdo et al.l 12013 ). 
A distance D ^ 1 kpc was estimated from an empir- 
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Figure 1. Field of J0633 as seen with Chandra/ACIS-S in the 
0.3-8 keV range. The image is binned by 4 ACIS pixels and 
smoothed with a 4-pixel Gaussian kernel. The pulsar is marked 
in the image and its PWN adjacent to the pulsar from south is 
clearly seen. 60" x 70" dashed box shows the region used to ex¬ 
tract the PWN spectrum. The unrelated background source ‘A’ 
falling in the extraction region is also marked. The solid rectan¬ 
gle with dimensions of 60" x 70" shows the region used for the 
background extraction. The intensity is given in counts per pixel. 


i cal relation between 7 -ray and spin-down luminosities 
( Saz Parkinson et al ] 12mft . In addition, J0633 powers 
an extended pulsar wind nebula (PWN) visible in X- 
rays south fr om the pul s ar. A nalyzing 20 -ks Chandra 
observations, iRav et aP ( 2Qll[) found that the X-ray 
spectrum of J0633 contains a thermal component which 
dominates at low energies and a non-thermal power-law 
(PL) component of the NS magnetospheric origin de¬ 
scribing the high energy spectral tail. 

Re-analysing the Chandra data, we found a hint of 
an absorption feature in the spectral fit residuals. We 
argue that this feature is real in Section 12.11 and dis¬ 
cuss its possible nature in Section 13.11 In Section 12.21 
we analyse the X-ray spectr al continuum. We confirm 
findings of iRav et al l ( 2011 ) and extend their analysis 
by incorporating natural constraints on the interstellar 
absorption and the distance to the pulsar. We analyse in 
detail the thermal component and investigate whether 
it can be attributed to the emission from the entire NS 
surface or a substantial part of it. In such a case, it is 
possible to confront inferred temperatures with predic¬ 
tions of NS cooling theories; this is done in Section [321 
In addition, a speculative birthplace of J0633 is pro¬ 
posed in Section 13.41 If it is real, it provides additional 
independent constraints on the distance to the pulsar. 


2 ANALYSIS OF THE X-RAY DATA 

We retrieved the dat 80 from the Chandra archive. Data 
mode was VFAINT, exposure mode was TE and the 
pulsar was exposed on the ACIS-S3 chip. The CIAO 
v.4.6 chandrajrepro tool with CALDB v.4.5.9 was 
used to reprocess the data set. A fragment of the Chan¬ 
dra image of the pulsar field is shown in Figure [H The 
pulsar and its extended PWN are clearly seen in the im¬ 
age. We extracted spectra of the pulsar and the PWN 
in the range of 0.3-10 keV with the CIAO v.4.6 specex- 
tract tool. For the background, we used a region free 
from any sources which is shown by the solid rectangle 
in Figure [T] The PWN spectrum was extracted from 
the dashed rectangle shown in Figure [T] excluding the 
pulsar and the point-like background object ‘A’ which 
overlap with the PWN. The number of counts for the 
PWN and the background in the same region were 397 
and 402, respectively. To extract the pulsar spectrum, 
we used a circular aperture centred at the pulsar po¬ 
sition with the radius of 2 '.'5, which ensures maximal 
signal-to-noise ratio. There are 332 pulsar counts (> 
98% of the total number of pulsar counts), two counts 
of background and two counts of the PWN within the 
aperture. 

We fitted the pulsar spectrum by an absorbed sum 
of power-law (PL) and therr nal cornponen ts using 
the XSPEC v. 12.8.2 package ( ArnaudI 1 19961 ). To ac¬ 


count for the photoelectric absorption we used the 
XSPEC PHABS model with default cross- sections bcmc 
( Balucinska-Church fc McCammonI 19921 ) and abun¬ 
dances angr (| Anders fc Grevess^ Il989[ ). For the ther¬ 
mal component we tried blackbody (BB) and hydrogen 
magnetic atm osphere niodels NS A (jPavlov et al.lll995[ ) 
and NSMAX ( Ho et al.1l2QQ8[) . Since the number of the 
pulsar counts is small, binning the spectrum goes at the 
expense of spectral resolution. Therefore, we used un¬ 
binned pulsar spectrum i n our analys is. Accordingly, we 
employed the C-statistic (|Cashlll979l ) for fitting, instead 
of more common 

We performed the fitting by a Markov chain Monte- 
Carlo (MCMC) sampling procedure assuming uni¬ 
form prior distribution for model parameters. We 
employe d the affine-invarian t MCM C sampler devel¬ 
oped by [Goodman fc Wear3 ( 2010 ) and implemented 
i n a p ython package emcee bv [Foreman- Mackev et ar 
( 2013 ). For each model we used a set of 100 MCMC 


walkers performing 1500 steps after initial burning, 
which is large enough c onsidering that typical an to- 
correlation time (see, e.g.. [Goodman fc Wearel[ 201 (i for 
details) was of the order of several tens (50-80). This 
resulted in a set of 150000 samples in total which was 
enough to reliably approximate the posterior distribu¬ 
tion of the model parameters. Having the sampled pos¬ 
terior distribution we obtained best-fit estimates and 

^PI Roberts, Chandra/AClS-S, Exp. time 20 ks, OBsID 11123 
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credible interval^ of the model parameters, and corre¬ 
sponding values of C-statistic. 

To assess the goodness of fit for each model we simu¬ 
lated spectra under the model in question with param¬ 
eters drawn from the corresponding sampled posterior 
distribution. Fitting these spectra we obtained the ref¬ 
erence distribution of the C-statisticH The C-statistic 
value obtained for the real data was compared with this 
distribution. Goodness-of-fit tests performed in such a 
way showed that any one-component model, pure PL 
or pure thermal, fails to describe the data (for 100% of 
simulated spectra, the C-statistic is less than observed 
one). In contrast, all two-component models mentioned 
above fit the data well, that is, the observed C-statistic 
value in each case is within one standard deviation from 
the mean of the corresponding reference distribution. 
Similarly, we found that an absorbed PL model is con¬ 
sistent with the spectrum of PWN. 


2.1 Absorption feature 

The pulsar spectrum is shown in the top-left panel of 
Figure [2] together with the best-fit PHABSx (BB-fPL) 
model. The C-statistic value of 266.9 is also presented 
in the plot. The corresponding fit residuals are shown in 
the bottom left panel of Figure [2l While the unbinned 
spectrum is used for fitting, the hard-energy part of 
the spectrum in Figure [2] is binned for illustration pur- 
posesEI 

What attracted our attention, was a hint of an ab¬ 
sorption feature in the fit residuals at about 0.8 keV. 
The approximate feature position is marked by the thick 
bars in the top panels of Figure O There are at least five 
consecutive channels which seem to stay apart from the 
best-fit model. We thus added the Gaussian absorption 
(GABS) component to the model and refitted the data. 
The pulsar spectrum with the new best-fit model and 
corresponding fit residuals are shown in the top-right 
and bottom-right panels of Figure [21 respectively. The 
new fit gives better C-statistic value of 254.1. Similar 
effect was observed for all continuum models we had 
testedH 


The difference in the C-statistic values between the 
models with and without the line is AC = 266.9 — 
254.1 = 12.8. In order to estimate the statistical signifi¬ 
cance of the fit improvement, we constructed the appro¬ 
priate reference distribution for AC, or likelihood ratio 
test (LRT) statistic, in a manner similar to the pro¬ 
cedure of assessing the goodness of fit. We simulated 
spectra under the model without line (the null model), 
drawing parameters from the corresponding posterior 
distribution sampled via MGMG. We then fitted simu¬ 
lated spectra with the null model and the model with 
line (the trial model) and computed the corresponding 
AC or LRT statistic. The LRT distribution basing on 
5000 data sets simulated with BB+PL as a null model 
is shown in Figure |3l where AC = 12.8 obtained for the 
data is shown by the vertical dashed line. It is seen, that 
only for 9 out of 5000 simulations the improvement in 
C-statistic was greater than 12.8. This means that such 
an improvement can hardly happen by chance if the null 
model is the true one. This statement can be quantified 
by the posterior predictive p-value, that is, a fraction of 
simulations with the LRT larger than the observed one. 
In our case p-value is 0.002 which favours the absorption 
line presence. Similar analysis performed with other 
continuum models resulted in p-values of the same or¬ 
der of magnitude. This method is known as a method of 
posterior predictive p-values, th e Bayesian mode l check - 
ing approach r ecommended by jProtassov et al.l (|2002[ ) . 
In particular, jProtassov et al. ( 20021) argue that this 
method is superior to the more common F-test in as¬ 
sessing the presence of additional model coni ponent. We 
refer the reader to, e.g.. lGelman et al.l (|2003[ ). for a text¬ 
book description of the method. 


Table 1 Median values of the absorption feature parameters with 
BB+PL as a continuum model. 


Eo (eV) 

a (eV) 

T (eV) 

EW (eV) 

8041^2 

<285 

>10 



90% credible intervals for the line centre Eq and equivalent width 
EW are given, while 99.9% limits for the line depth r and the 
width a are presented. 


^The credible interval is any continuous part of the parameter’s 
marginal distribution containi ng certain fraction of the total dis¬ 
tribution (|Gelman et al.l2003l) . Here we adopt the range between 
the 5% and 95% quantiles, unless stated otherwise. 

^This procedure is close to what the XSPEC goodness task does. 
The goodness task uses the standard bootstrapping scheme 
where the best-fit estimates of the model parameters are used to 
simulate spectra. The approach we employed is a more general 
as it incorporates the model parameter uncertainties conditional 
on the current observational data. 

"^We keep spectrum unbinned up to channel 64, group 8 channels 
into 1 bin for channels 64-128 and 16 channels into 1 bin for 
channels 128-1024. 

^This result would be the same with any model for continuum 
which is smooth across the putative line region. 


The best-fit spectral line parameters and their uncer¬ 
tainties extracted from MGMG are presented in Table [H 
Here we use the Gaussian line model which contains 
three parameters and is given by the following expres¬ 
sion 

GABS(E)=exp(-^e-^^), (1) 

where E is the photon energy, Eq is the line centre, a is 
related to the full width at half maximum (FWHM) of 
the line as FWHM ^ 2.35cr, and parameter r regulates 
the line depth. Then the optical depth at line centre 
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E[keV] E[keV] 

Figure 2. Chandra/ACIS-S spectrum of J0633 fitted by an absorbed BB+PL (left) and by an absorbed (BB + PL) x GABS models 
(right). Best-fit models are shown with black solid lines in the top panels. Red and cyan lines show BB and PL model components, 
respectively. The fit residuals in form of cire shown in the corresponding bottom panels. The absorption line position is shown by the 
thick bars. Best-fit C-statistic values are shown for both fitting models in the top panels with the values of the fit degrees of freedom 
given in parentheses. 



Figure 3. Probability density function (p.d.f.) for the likelihood 
ratio test (LRT) statistic, that is, the difference in the (7-statistic 
for the BB+PL and (BB + PL) x GABS fits for 5000 simulated 
data sets. Vertical dashed line indicates the observed LRT statis¬ 
tic A(7data = 12.8. The corresponding p-value is also shown, see 
text for details. 

is a more direct measure of a line strength 

is the equivalent width (EW) which is defined by the 
following expression 

+ 00 

EW = J (1 - GABS(E)) dE (2) 

— OO 

The main advantage of EW is that it is weakly depen¬ 
dent on the particular shape of the spectral feature. Eor 
the Gaussian line in the optically thin regime (r/cr ^ 
1), EW r. 


In Eiguredl we show one- and two-dimensional (ID 
and 2D) marginal posterior distributions for the line 
parameters. As seen, Eq is well constrained around 0.8 
keV. The median value of Eq along with the 90% cred¬ 
ible interval is presented in Table [H Unfortunately, the 
situation is different for the line width and depth. Eig- 
ure m reveals a bimodal, worm-like, 2D posterior distri¬ 
bution for a and r parameters. This bi-modality is also 
seen in the ID distributions for these parameters. The 
“worm” head and body correspond to different types of 
the absorption line. The worm-body mode corresponds 
to a strong saturated line with the width smaller than 
the Chandra/AClS-S spectral resolution (EWHM ^ 100 
eVEI) . The fit quality is then determined by the wings of 
the line, which results, as can be shown analytically, in 
a strong degeneracy between the line width and depth 
giving the long worm-body valley in the likelihood dis¬ 
tribution. On the other hand, the worm head corre¬ 
sponds to a broader and weaker line. With the present 
data we cannot discriminate between these possibilities. 
Therefore, only the upper limit on a and lower limit on 
r are given in Table [TJ At the same time, EW is well 
constrained as seen from Eigure |4] and Table [T] which 
can be regarded as the most straightforward argument 
in favour of the line. 

Note that the above results are almost independent 
on the particular continuum model used to fit the pul¬ 
sar spectrum. In addition, they remain qualitatively the 
same if models other than GABS are used to fit the ab¬ 
sorption feature, for instance the models for a cyclotron 
absorption line (CYCLABS in XSPEC) or an ionization 
edge (EDGE in XSPEC). 
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Figure 4. ID and 2D marginal posterior probability distributions for the line parameters (line central energy Eq^ width cr, depth r 
and equivalent width EW) in the (BB + PL) x GABS model. The 5%, 50% and 95% quantiles are shown with vertical dashed lines in 
the ID distributions. 


The Bayesian analysis shows that the chances are 
low that the absorption feature is caused by Poisson 
fluctuations of the data counts. It may be an instru¬ 
mental artifact though. It would thereby be appropri¬ 
ate to examine if there are similar features in spectra 
of other sources in the Chandra/ACIS-S field of view. 
However, all other point sources are substantially dim¬ 
mer than the pulsar, showing no more than several tens 
of counts, and analysis of their spectra is not conclu¬ 
sive. We also examined the spectrum of the PWN (see 
Figure [5]). Unfortunately, the PWN spectrum is more 
noisy at the soft energies than the pulsar spectrum due 
to much higher background. There is no line seen here, 
at least at the first sight. Indeed, we got improvement in 
statistic of only AC ^1.5 after fitting the PWN spectra 
with the PL x GABS model with Eq of about 0.8 keV. 
The posterior-predictive p-value test gives no evidence 
against the featureless continuum model with p-value ~ 
0.43. We also checked that there were no flares during 
observations which could distort the spectrum. Accord¬ 
ingly, the background spectrum is in agreement with 
the quiescent spectrum of the diffuse soft X-ray back¬ 


ground as seen with Chandra/ACIS (|Markevitch et al 


20031). 


2.2 X-ray continuum 


In this section we employ the same Bayesian technique 
to analyse the properties of the X-ray continuum as¬ 
suming the presence of the line. We will not explicitly 
indicate this further for brevity. 

Main problems arising during the analysis of the X- 
ray emission from NSs come from unknown distances 
D and interstellar absorption towards these objects. In 
order to get better constrains on the latter, we now fit 
simultaneously the pulsar and the PWN spectra, tying 
the value of the hydrogen column density A^h between 
the fits. Recall that the PWN spectrum is well described 
by an absorbed PL, while the pulsar spectrum contains 
thermal and non-thermal components. For the thermal 
component, as already mentioned above, we tested sim¬ 
ple blackbody as well as several magnetised hydrogen 
atmosphere models. All latter give similar results, there¬ 
fore we selected a par ticular model fr om NSMAX fam¬ 
ily, labeled 1260 (see IHo et al. 2008 ). The reason for 
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Figure 5. Chandra/ACISS spectrum of the J0633 PWN fitted 
by an absorbed PL. The spectrum is binned to ensure > 1 counts 
per bin. The absorption line position is shown by the thick bar. 


selection of this model is twofold. First, NSMAX mod¬ 
els account for the partial ionization in the NS atmo¬ 
spheres which makes their usage more physically mo¬ 
tivated for low temperatures in comparison with older 
NSA models. Second, the 1260 model corresponds to 
the surface magnetic field of 5 x 10^^ G which is close 
to the J0633 value as inferred from P-P observations 
assuming dipole losses. Any atmospheric model depends 
on the NS surface gravity. In the NSMAX models, this 
is incorporated via the gravitational redshift parameter 
1 -h 2 :. In our fit we fixed 1 -h 2 : = 1.21 which corresponds 
to a reasonable NS model with a mass Mns = 1.4M0 
and a circumferential radius i?NS = 13 km. The appar¬ 
ent NS radius as seen by a distant observer in that case 
is i? = (1 + z)R^s ~ lb km. We have checked that the 
redshift parameter is not constrained by the data and 
does not influence the final results. We also note that, 
due to effects of general relativity, the NSMAX model 
can be applied only for describing the emission coming 
from the entire surface of the NS, while the blackbody 
model can be used to describe the emission from any 
part of the stellar surface. 

In first two rows of Table [21 we show best-fit param¬ 
eters of the BB+PL and NSMAX-fPL models with un¬ 
certainties corresponding to 90% credible intervals. The 
latter were inferred from marginal Bayesian posterior 
distributions as descr ibed above . Thes e results are con¬ 
sistent with those of iRav et al ] (1201 ih and we thereby 
confirm their findings. The values of C-statistic are also 
shown in Table [21 

It is instructive to consider which part of the poste¬ 
rior parameters distribution can describe the physically 
allowed situation (in the selected model framework). In 
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particular, can the fit results for the thermal component 
correspond to the emission originating from the entire 
NS surface or a part of it? In Figure |6l we plot ID and 
2D marginal posterior distributions for the parameters 
of BB and NSMAX models, which are temperature 
T and the normalisation of the thermal component. For 
both models, the temperature is given as measured by a 
distant observer. The thermal component normalisation 
is presented in the form of the distance to the pulsar if 
the apparent radius of the emitting area is R = 16 km. 
As seen, the minimization of the likelihood leads to cor¬ 
relation between A^h and normalisation of the thermal 
component and thus to anti-correlation between A^h and 
distance D. In contrast, these latter quantities obviously 
correlate in nature. It is unlikely to have a high A^h value 
at a small distance and vice versa. 

There exist empirical models which provide depen¬ 
dence of interstellar absorption on distance. Such mod¬ 
els usually describe the correlation between distance 
and optical extinction Ay. The latter can be trans¬ 
formed to A^h using one of the empirical relations be- 
tween Ay and A/'h, for instance, the one given by 


Predehl fc Schmitt! (|l995f ). To obtain the relation be¬ 


tween Ay and distance in the pulsar direction, we made 
use of the three-dimensional m odel of Galactic extinc¬ 
tion from Drimmel et al. (l2QQ3l) and also took into ac¬ 


count the value of maximal Ny in the pulsar direc¬ 
tion, ^ (6-7) X 10^ ^ cm“^ ( Dickey fc Lockm^ll99Ql: 
Kalberla et al. I2QQ5I) . This information can be roughly 
summarized in a simple relation A^h [10^^ cm“^] « D 
[kpc] at D < 7 kpc. We present this relation in the D- 
A^h plate in Figure [6] with solid, dashed and dot-dashed 
lines assuming R = 16, I and 20 km, respectively. The 
latter value is a reasonable theoretical maximum for 
the NS apparent radius. Consider, for example, the BB 
model. As can be seen, 16 km radius is consistent with 
the data. The corresponding D-Nn relation crosses the 
marginal posterior distribution not far from its max¬ 
imum. The star in this case appears to be at about 
2.5 kpc from the Sun and has the temperature of about 
100 eV. If the apparent radius is I km, our analysis 
leads (dashed line) to A^h ~ 10^^ cm“^ and hence places 
J0633 at D « I kpc. In this case, the inferred temper¬ 
ature should be higher, more than 120 eV. The radii 
much lower than I km would give worse fits and are 
unlikely. The portion of the posterior distribution of 
parameters lying down-right from the dot-dashed line 
in Figure [6] corresponds to unrealistic R > 20 km. This 
means that although the fit is formally good there, the 
BB model with such parameters can not describe ther¬ 
mal emission from the NS surface. Looking at the N^-T 
and D-T plots we conclude that too low temperatures 
are impossible, however the corresponding regions are 
broad and this restriction is rather weak. Similar anal¬ 
ysis can be performed for the NSMAX model. In this 
case, most of the posterior distribution corresponds to 
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Table 2 Best-fit spectral parameters for continuum models. 


Model 

Nu 

(10^^ cm-2) 

r 

4- psr 

Rpsr 
( 10 “® ph 
keV“^ cm~^ 

T 

(eV) 

R “ 
(km) 

D “ 
(kpc) 

r 

-L pwn 

Rp'wn 
( 10 -® ph 

keV“^ cm~^ s~^) 

C/d.o.f. ^ 





No prior 






BB+PL 

2-411:4 

1 0 + 0-6 

9-611:? 

105111 

2Atl\t 

7+12 
‘ -5 

-i 9+0.3 

-^•^-0.3 

26.717^5^ 

381.7/792 

NSMAX+PL 

2.911:5 

1 4 + 0-6 

6-7I3.5 

4lll? 

0^ + 205 

0 i^3+^ 08 

U.OO_Q 45 

1 3+0.4 

-^•'^-0.3 

29.61^^4® 

388.3/792 





With prior 





BB+PL 

2-211:5 

1 0 + 0-6 

9 3+6.6 

1081^4 

5ir 

9 -| +2.2 

-i 9+0.3 
-^•^-0.3 

26.ll®J 

383.2/791 

NSMAX+PL 

1 7+0.6 
‘ - 0.7 

-1 9 + 0.6 

4 

^^ — 2.6 

53 ^ 7 ^ 

1215 

-^•'^-0.6 

-^•-^-0.2 

23.31?:? 

4.04.AI791 


Temperatures T and emitting area radii R are given as measured by a distant observer. Redshift parameter for NSMAX models is fixed 
at 1.21. r and K are the photon index and the normalisation of the PL component. All errors correspond to 90% credible intervals 
derived via MCMC. For models in two last lines an informative prior which includes information on the N^-D correlation is applied, 
see text for details. 

“ In the “no prior” case, R is given assuming D = 1 kpc and D is given assuming R = IQ km. 

^ The number of degrees of freedom is different from those given in Figure [2] because here the PWN spectrum is included. 


the region with i? > 20 km. However, radii of the order 
of 10 km are still allowed, giving A^h ~ 1.5 x 10^^ cm“^ 
and D ^ l.b kpc. This shrinks the possible tempera¬ 
ture range (in contrast to BB model). As clear from 
Figure [6l temperature is constrained at T > 40 eV and 
A^h < 2 X 10^^ cm Again, the low-temperature part 
of the posterior distribution requires too high emitting 
area radii. Finally, as seen from the plot, the 1 km ra¬ 
dius is too small if the NSMAX model is applied. 

The Bayesian approach provides natural framework 
for inclusion of the additional information, such as the 
D-Nn relation discussed above, by defining the appro¬ 
priate prior distribution. Moreover, it easily allows to 
take into account uncertainties in the prior knowledge. 
Using this possibility, we incorporated the D-N^ re¬ 
lation as the Bayesian prior in the following way. We 
made a conservative assumption that the relation A^h 
[10^^ cm“^] D [kpc] is accurate up to a factor of 
two. In principle, the central values are more likely, and 
a bell-shaped form of the prior, for instance, Gaussian 
would be appropriate. However, as the exact variances 
of the D-Nii relations are unavailable for us, we used 
the following flat prior: O.bD [kpc] < A^h [10^^ cm“^] 
< 2D [kpc]. In addition, we constrained D to be less 
than 7 kpc, which is the approximate distance to the 
edge of the Galactic disk in the pulsar direction. The 
application of such prior naturally splits the thermal 
normalisation into two independent parameters, D and 
R. We also constrained the latter to be less than 20 km. 

The best-fit parameters inferred with account for the 
prior are shown in third and fourth rows of Table [2l 
Both BB and NSMAX models pass the goodness-of-fit 
test. The marginal ID and 2D posterior distributions 
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for A^h, T, D and R are shown in Figure [71 We can see 
how the prior works, comparing Figure [7] with Figure |6l 
The distance is now constrained by the prior relation 
and A^h, and the constrains on R are obtained from 
D and the thermal component normalisation. The pa¬ 
rameters in Table [2] generally agree with the qualitative 
considerations presented above. We see that both BB 
and NSMAX models are consistent with the physical 
picture where emission comes from the entire surface of 
the star. For the BB model, however, the star should 
be somewhat more distant and more absorbed than in 
case of the NSMAX model. The hot-spot (of about 1 km 
size) interpretation is also not excluded under the BB 
model. As it usually happens, the best-fit BB temper¬ 
ature is about twice hi gher than that for a hydrogen 
atmosphere model fe.g.. jPavlov et al.ll2QQl[) . 


3 DISCUSSION 

3.1 Absorption feature 

The analysis performed in Section [2T] favours presence 
of the absorption feature in J0633 spectra at about 
0.8 keV. Unfortunately, the shape of the feature is 
poorly constrained with the current data, precluding us 
from plausible interpretation of its nature. The first pos¬ 
sibility to consider is the cyclotron line. The cyclotron 
absorption line position, as seen by a distant observer, 
for a particle of charge Z and mass m is given by 

TTi 

& = 11.577(1 + z)-^Z—Bi 2 keV, (3) 

m 

where rrie is the electron mass and the line is assumed to 
form at the NS surface. Now we can estimate a surface 
magnetic field as H 8 x 10^^ G if the line is produced 
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by electrons and B ^ 1.4 x 10^^ G if it is produced 
by protons, and even higher values for more massive 
ions. Both values are inconsistent with the spin-down 
estimate of the dipole magnetic field, B = 4.9 x 10^2 G. 
The “cyclotron” field is much lower in case of electrons 
and much higher in case of protons. 

Note that for other INS showing absorption lines, 
spin-down magnetic fields, when determined, usu¬ 
ally disagree with “cyclotron” magnetic fields. For 
SGR 0418+5729 and PSR J1740+1000, the discrep¬ 
ancy is as strong as for J0633. In these cases, 
the proton cyclotron line interpretation is possible, 
for instance, if there are strong small-scale (mul¬ 
tipolar) surface components of the magnetic field. 
The presence of the s mall-scale fields is widely dis¬ 
cussed in literature (e . g., Xsseo fc Khechinashvilill2QQ2[ 
Harding fc MnslimovI l2Qll[ and references therein). 
Tiengo et al. ( 2Q13h suggested this interpretation for 


the spectral feature of SGR 0418+5729. They, how¬ 
ever, had additional arguments from the phase-resolved 
spectral analysis which is unavailable in case of J0633. 
A somewhat similar feature was recen tly detected in 
the sp ectrum of INS RX J0720.4—3125 ( Borghese et al.l 


20151 ). The feature is at ^750 eV which corresponds, if 


interpreted as the proton cyclotron line, to the magnetic 
field about seven tim es higher than the spin- down one. 

On the other hand. [l^r gaitsev et al. ( 20121) proposed 
that the absorption feature in PSR J1740+1000 is the 
electron cyclotron line which is produced by a popula¬ 
tion of warm electrons occupying some regions in the 
pulsar magnetosphere similar to Van Allen belts in the 
Earth magnetosphere. If we adopt this interpretation 
for the line in J0633 and assume that the magnetic field 
approximately obeys the dipole law, B oc we can 
estimate the position of the magnetospheric radiative 
belt diS r ^ 4Rns5 oi* about 30-40 km above the NS sur¬ 
face. 

For the remaining objects, the disagreement is less 
dramatic. The GGO 1E1207—5209 shows at least two 
features in the X-ray spectrum. Its low spin-down mag¬ 
netic field suggests the electron cyclotron interpreta¬ 
tion. However, the position of the fundamental har¬ 
monic estimated from the spin-down value is larger by a 
factor o f 1.4 than the positio n of the strongest spectral 
feature (|Gotthelf et al.llioi^ . Finally, for INSs from the 
magnificent seven group, spin-down magnetic fields are 
lower by a factor of 1.1-7.2 than proton cyclotron mag¬ 
netic fields estimates, with the best agre ement achieved 
for RX J1308.6+2127 (IPires et al.ll2014h . 


Another possible explanation of the absorption fea¬ 
ture suggests that it results from atomic transitions 
either in the stellar atmosphere or in the interstellar 
medium. Varying abundances in the model for the inter¬ 
stellar photoelectric absorption we found that the fea¬ 
ture could be explained assuming overabundance of Fe 


along the pulsar line of sight. It is, in principle, pos¬ 
sible, as the J0633 position projects onto the Mono- 
ceros Loop nebulosity, see Figure [H which was rec- 


tical, and X-rays ( 

IDaviesI 

1963 

[ iGebel & Shore! 1972: 

iDavies et al. 19781: 1 

Leahv et al. 1 

49851). The distance to 


the remnant is not exactly known, however most of 
the e stimates suggests a value of approx imately 1.6 kpc 
(e.g.. lBorka Jovanovic fc Urosevidl2QQ9[ and references 
therein) providing an SNR shell diameter of about 
0.1 kpc. This distance together with our estimates for 
the J0633 distance (Table [2]) allows the pulsar to be 
behind the SNR and suffer from an additional absorp¬ 
tion of the SNR origin. However, to obtain a good fit 
for the J0633 spectrum, too large Fe/H ^ 3 x 10“^ is 
required, which is about ten times as large as the solar 
one. To provide such a high abundance of Fe on average 
along the pulsar line of sight, abundance of Fe in the 
Monoceros Loop SNR itself should be at least 100 times 
the solar abundance, which seems implausible. Addi¬ 
tionally, such a strong modification of the interstellar 
absorption should also affect the PWN spectrum, how¬ 
ever, as stated above, we did not find any signature of 
spectral features there. Finally, if the absorption feature 
is formed in the mid-Z atmosphere of NS, th en generally 
broad er and weaker features are expected (|Mori fc Hoi 


20071) . 



209 208 207 206 205 204 203 202 201 200 

Galactic longitude 


Figure 8. Ha image of the Monoceros Loop region in Galactic 
coordinates t aken from the South ern H-Alpha Sky Survey At¬ 
las: H-Alpha (|Gaustad et al.|[200lh . The Monoceros Loop SNR is 
marked by the dashed circle. The position of the pulsar and its 
possible proper motion direction are shown by X and the arrow, 
respectively. The Rosette nebula suggested as its likely birthplace 
is pointed on. 
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Figure 9. Observations of isolated cooling NSs vs. cooling the¬ 
ory predictions. Temperatures obtained utilising the BB model 
are shown with the cyan colour, and those obtained with various 
atmospheric models are shown with the blue colour. The J0633 
data points for BB and NSMAX models are shown by the star 
symbols. Dashed lines present the cooling curves corresponding 
to standard cooling of 1-1. 9Mq NSs with the APR EOS. The 
filled region corresponds to a possible range of standard cooling 
curves including the unrealistically compact equations of state. 
Upper and lower solid curves illustrate the effects of the nuclear 
superffuidity in the NS core. Upper curve correspond to strong 
proton superffuidity suppression of the cooling, while lower curve 
is calculated including also Cooper pair formation emission from 
the triplet neutron superffuid. See text for details. 


3.2 J0633 in the view of the NS cooling 
theories 

According to the Tableland Figure [3 the J0633 ther¬ 
mal emission can originate from the entire stellar sur¬ 
face. In this case, it is instructive to compare the re¬ 
sults with the NS cooling theories. It is worth doing, 
according to Table [21 for both BB and NSMAX models, 
although the inferred surface temperatures are differ¬ 
ent. This comparison is performed in Figure |9l where 
the positions of J0633 on T-r plane are shown for 
both models along with the data for other cooling iso¬ 
lated NSs. The data on the la tter objects are taken 
from the references c o llected by [Yakovlev et ah! (2008) 


and Kaminker et al 


, excluding upper limits, 
with addition of several sou rces. The addit i onal s ources 
include PSR .71741 -2054 (iKarnova et al.l Eflill . PSR 
J0357+3205 (iKiricheuko et a].ll2ni4. 1^ J0007+7303 
in the CTA 1 SNR (|Ca,ra.veo et alJbni fi ITju et a,l1l2ni flh 
and two CCOs, where carbon atmosphere is used to 
describe their th ermal radiation. Thes e are CCO in 
the Cas A SNR ( Yakovlev et al.l [201l[) . the youngest 
source in Fignrepand CCO XMMU J173203.3-344518 


(Klochkov et al 


2015[ ). the hottest source in Figure [9l 


For the latter source we combined temper atures for 
both distances given in iKlochkov et al.l(|2015h . and mul¬ 
tiplied the temperature errorbars by a factor of 2 to es¬ 
timate 2(7 errors. In Figure [H with cyan colour we show 
point obtained with BB model, and with blue colour 
points obtained with various atmospheric models, see 
the above cited references for details. In Figure [9l we 
also artificially adopt a factor of two uncertainty on 
J0633 age. It is seen, that if J0633 is covered by the 
atmosphere, it is the one of the coldest middle-aged 
(r < 10^ yr) isolated NS with measured surface tem¬ 
perature. If, in contrast, the BB model actually fit the 
data, the inferred temperature is much higher. 

According to theory, isolated NSs cool down via 
the neutrino emission from their interiors and via the 
photon emission from their surfaces. The middle-aged 
stars are of particular interest as they have isother¬ 
mal interiors, except a thin heat blanketing layer near 
the surface, and the ir cooling is dominated b y neu¬ 
trino emission (e.g., Yakovlev fc PethickI l2QQ^ . Mea¬ 
surements of surface temperatures of such stars allow 
one to determine the neutrino cooling rate and there¬ 
fore to directly explore the properties of matter deep 
inside the star. By the filled region in Figure [9] we show 
the predictions from the so-called standard NS cooling 
scenario which assumes that a star has nucleon core 
and cools mainly via the mod ified Urea processes o f 
neutrino emission. According to [Yakovlev et al.[ (|2Qllh . 
the standard cooling mainly depends on the compact¬ 
ness of the star x = Rg/R^s^ where Rg is the gravita¬ 
tional radius, or equivalently on the gravitational red- 
shift 1 z = (1 — being largely independent on 

the particular NS model. More compact stars gener¬ 
ally cool faster. This property allows to estimate the 
neutrino cooling rate of a particular middle-aged star 
without performing cooling simulations. The filled re¬ 
gion corresponds to a broadest possible region that can 
be reached with the standard cooling (note that here 
the standard iron heat-blanketing envelope is used). 
It includes also the cooling curves corresponding to 
highly unrealistic equations of state, which allow for 
extre mely compact stars (for details see [Yakovlev et af 


l2Qll[) . For comparison, in Figure [9] with short-dashed 
lines we present the cooling curves for stars with a par- 


tion (same as used by, e.f 

1 .. Yakovlev et al. 20111 of the 

APR EOS (Akmal et al. 

|l998l 

), widely used as a stan- 


dard. The curves are given for a range of NS masses 
from 1.0 to 1.9M0 per 0.1 solar mass, plus the curve 
for the maximal mass Mmax = 1.929M0 for this EOS. 
The direct Urea process is, in principle, allowed in the 
massive stars with APR EOS, but here it is switched 
off. The standard cooling curves for other reasonable 
EOSs basically fall in the same region. In other words, 
the part of the filled region which is colder than the 
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lowest dashed curve is reached in the standard cooling 
scenario in principle^ but is marginal. 

Foll owing the method described in I Yakovlev et al.l 
dmH), we found that the neutrino cooling rate of J0633 
should be 30-1000 times stronger than the standard 
one, if the NSMAX model is applied, for a reasonable 
star compactness x < 0.5, and with account for a factor 
of two uncertainty in the pulsar age. Only for unrealisti¬ 
cally compact stars, x ~ 0.7, the inferred NSMAX tem¬ 
perature can be reached in the standard cooling mod¬ 
els. A moderate increase in the neutrino cooling rate 
<100 can be explained by the minimal cooling theory 
which includes also the neutrino emission in the pro¬ 
cess of Cooper pair fo rmation in triplet neutron super¬ 


fluid in the NS core ( Page et ^l2004 ICusakov et al 


20041 ). We illustrate this possibility with lower thin 


solid line in Figure [9l It corresponds to a I.TMq APR 
EOS star and simila r superfluidity model as used by 
Shternin et al ] (l201lh in explanation of the data on the 


NS in Cas A SNR. Too lowering of temperature is hardly 
possible in the minimal cooling scenario. Nevertheless, 
in any case low temperatures of cooling NSs can be 
explained if the direct Urea processes are allowed in 
their cores. However, in order to get the temperatures 
like J0633 has, assuming NSMAX spectral model, these 
processe s have to be suppressed, fo r instance, by super¬ 
fluidity (|Yakovlev fc Pethi^ 20041 ). Otherwise the en¬ 
hancement of the neutrino emission will be too strong. 

For the BB model, the neutrino cooling rate, in con¬ 
trast, should be much weaker. We find that it must be 
suppressed by a factor of 10-300. The lowering by a 
factor < 50 is also possible in the minimal cooling sce¬ 
nario if strong proton superfluidity is involved which 
suppresses the conventional mUrca processes, and if the 
internal temperature of the star is hotter than the neu¬ 
tron superfluidity critical tempe rature so that Coope r 
pair emission does not operate (|Gusakov et al.N2004l) . 
This is illustrated in Figure [9] with the upper solid line 
which corresponds to IMq star with the APR EOS 
and strong proton superfluidity in the core. This curve 
fits nicely the BB data. Another possibility allowed in 
the cooling theory is that the heat-blanketing envelope 
of the star contains sufficient amount of the light ele¬ 
ments. Then the envelope is more transparent to heat 
and the star looks hotter than the star with the same 
i nternal thermal state , but iron (non-accreted) envelope 
( Chabrier et al.lfToQTl) . However, in the latter case the 
hydrogen (or other light-element, for instance, carbon) 
atmosphere would be more appropriate than BB to de¬ 
scribe the emission spectra. Also, the star will look hot¬ 
ter if magnetic field as strong as > 10^^ G is prese nt in 
the heat-blanketing envelope (jPotekhin et al.ll2QQ3[ ). Ei- 
nally, some additional heat ing mechanisms can opera te 
in the stellar interiors fe.g. I Yakovlev fc PethickI 2004 ). 

The inferred parameters of the BB model allow a dif¬ 
ferent interpretation of the NS thermal emission. It is 


possible that it actually comes from a hot spot on a 
colder NS surface which is heated by charged parti¬ 
cles coming from the magnetosphere along the mag¬ 
netic field lines near the magnetic poles. The conven¬ 
tional polar cap radius for J0633 can be estimated as 
^cap = 0.145 (Rns[ 10^ km])^/^(P[s])“^/^ km « 400 m. 
This is inconsistent with the spectral fit results (Ta¬ 
ble [2]) , however the 1-2 km emitting area radii are pos¬ 
sible. This range of radii correspond, according to Eig- 
ure[71 to the temperatures > 125 eV and smallest possi¬ 
ble distances of about 1-1.5 kpc. Note that these values 
are favoured by the pseudo-distance relation, and also 
by the possible birthplace of the pulsar in the Rosette 
nebula, see below. In the hot spot picture, inferred tem¬ 
perature cannot be compared with the cooling theories. 
The bulk of the surface is then cold and invisible in 
X-rays. 


3.3 Non-thermal efficiencies and luminosities 


The distance ranges inferred from the spectral fits (Ta¬ 
ble [2]) allow to constrain the non-thermal luminosities 
of J0633 and its PWN. In Table[3]we give the 2-10 keV 
X-ray fluxes of the PWN, and the non-thermal 

(PL) spectral component of the pulsar, , for both 
spectral models usec0. As expected, these values almost 
do not depend on the type of the thermal continuum 
model (BB or NSMAX). Corresponding non-thermal 
luminosities and are also given in Table [3] 

along with the values of efficiencies /E and 


^pwn _ These values show some dependence 


on a choice of the spectral model because in the two 
models the inferred distance ranges are slightly differ¬ 
ent (Table [2]). In any case, the parameters of the X- 
ray non-thermal emission in Table [3] are not peculiar 
in comparison with those for o t her pulsars with simila r 
E ( Kargaltsev fc PavlovI I 2 QQ 8 : Danilenko et al.l [20131) . 
In addition, in Tabled we show J0633 y-ray luminosi¬ 
ties and corresponding efficiencies 77 ^®^ = /E. 

They are calculated basing on the pulsa r y-ray flux 
(9.4 ± 0.5) X 10“^^ erg cm“^ s“^ ( Abdo et~al] 
I2Q13I ). Note, that for large distances, D > 3 kpc, al¬ 
lowed for the BB+PL model (Table [3j), the y-ray effi¬ 
ciency is higher than 1. However, there are y-ray pulsars 
with precise l y me asured distances which have 77 ^ > 1 
(|Abdo et al.ll2Q13l) . The apparent violation of the en¬ 
ergy conservation law is usually resolved by account for 
unknown beaming of the y-ray emission. 


3.4 Presumed birthplace 

The J0633 PWN morphology and extent (Eigure [T]) is 
reminiscent of, for instance, a well-studied bowshock 


^Two models with prior from two last rows in Table [2] 
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Table 3 Non-thermal fluxes, luminosities and efficiencies. 


Model 

logFP- 
(erg cm ^ s 

log 

(erg s“L 

log ryP®"' 

logFP-" 

(erg cm“^ 

logLP-" 
(erg s“^) 


logLP- 
(erg s“^) 

log<^ 

BB+PL 



o ' 7 + 0.6 
‘ - 0.8 


09 1 +0.6 

— 3 O'*”®'® 
'^•'^-0.8 

34.77°;« 

-O'+o'.s 

NSMAX+PL 

-13.31“;^ 

3i.oi“:= 

_4 .+ 0.5 


31.7l“J 

_3 4+0.5 

'^•^-0.5 

34.3t°-j 



X-ray fluxes, luminosities and efficiencies are calculated in the 2-10 keV range. 


“Mouse” nebula (G359.23—082 ) powered by the fast- 
moving PSR J1747-2958 (e.g., |HaieseE^|2^. The 
similarity with the Mouse suggests the direction of a 
proper motion of J0633 as shown by the long arrow 
in Figure [TJ Adopting a typical synchrotron cooling 
time of X-ray emitting electron s in PWNe of ^ 1000 
yr (e.g., iKargaltsev et al.l l2008l ) and the J0633 PWN 
tail size of 1.'3 (see Figure [1]), we estimated the pul¬ 
sar proper motion as 80 mas yr“^, suggesting the pul¬ 
sar angular displacement of 1?3 during its lifetime of 
^ 60 kyr. Looking backwards the assumed proper mo¬ 
tion direction on the extended field of view (Figure [8j), 
we found a likely birthplace of the pulsar at the esti¬ 
mated angular displacement - the Rosette nebula. It is 
known as a young, 50 Myr, active star forming region 
located at the edge of the Monoceros Loop SNR. The 
distance to the Rosette of ^ 1.5 kpc (|Ogura fc Ishida 


19811) is compatible with the J0633 distance estimates, 


discussed above. Looking from the other side, adopting 
the Rosette nebula as a likely birthplace we can inde¬ 
pendently estimate the distance to the pulsar. Taking 
into account the conservative uncertainties of a factor 
of two for the pulsar age and angular distance to its 
specific birth location inside the Rosette nebula, and as¬ 
suming that the pulsar 3D veloc ity is smaller than 2000 
km s“^ (e.g. Hobbs et al.l 1200^ . we obtain a distance 
range of 1.2 < D < 1.8 kpc adopting the distance to the 
Rosette nebula of 1.4-1.6 kpc. This range is consistent 
with that derived from the spectral analysis and does 
not put additional constraints on the thermal emission 
models. 


4 CONCLUSIONS 

We have analysed in detail the X-ray spectrum of the 
7 -ray pulsar J0633-f0632. We found the evidence of the 
narrow absorption feature in its spectrum at 804^26 
with the equivalent width of 63'^^q eV, where errors are 
at 90% confidence. While the shape of the feature can¬ 
not be constrained with the current data, the failure of 
any smooth continuum model in vicinity of 0.8 keV is 
statistically proven. We briefly discussed possible physi¬ 
cal interpretations of the detected feature and currently 
favoured the cyclotron nature. 
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Apart from the spectral feature, we investigated the 
properties of the X-ray contin uum. We confirmed the 
conclusion bv iRav et al.l ( 201l[ ) that the soft part of the 
J0633 spectrum is dominated by the thermal compo¬ 
nent, and the hard tail is described by the non-thermal 
PL. We went further and took into account the corre¬ 
lation between the distance to the pulsar and the in¬ 
terstellar absorption along the J0633 line of sight on 
the basis of the empirical distance-extinction maps. It 
was included in the form of the prior distribution for 
the model parameters. In addition, the PWN spectrum 
was fitted simultaneously with the pulsar spectrum, al¬ 
lowing to better constrain the Nu parameter. As a re¬ 
sult of this analysis, we found that the thermal emis¬ 
sion possibly originates from the entire surface of the 
star and its spectrum can be equally well described by 
the blackbody or the magnetised hydrogen atmosphere 
models. In the BB case, the hot spot origin of the ther¬ 
mal emission is also possible. The distance to the pulsar 
was constrained within 1-4 kpc range basing on the X- 
ray spectral analysis. This is especially important, as 
the dispersion measure distance is unavailable for the 
radio-silent J0633. 

Confronting the inferred temperatures with the data 
on other cooling NSs and predictions from NS cool¬ 
ing theories, we found that for the atmospheric spec¬ 
tral model, J0633 is one of the coldest middle-aged NS 
with measured temperature. In this case it must cool 
considerably faster than the standard cooling scenarios 
predict. In contrast, if the spectrum of the substantial 
part of the NS surface is blackbody, then J0633 is hotter 
than a standard cooling NS. 

In addition, we found a possible birth site of the pul¬ 
sar - the Rosette nebula. This, along with the shape 
of the PWN, suggests that J0633 can have prominent 
proper motion. At the moment, these findings do not 
impose additional constrains on the pulsar properties. 
When this paper was in preparation, deeper observa¬ 
tions of J0633 with XMM-Newton were approved for 
the AO-14 observational cycle. If these observations are 
performed, more elaborate consideration of the absorp¬ 
tion feature origin and nature of the continuum emission 
will be possible. The detection of X-ray pulsations will 
be especially useful to attribute the thermal emission to 
a small hot region or to the entire surface of the NS. The 
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study of the variation of the absorption feature with the 
rotational phase is also important to explore its origin 
and, if the feature is a cyc l otron line, magnetic field 
geometry (IKargaltsev et al.l I2Q12I: IXiengo et al.l 12013 : 


Borghese et alJ[2Ql^ . The observations will also allow 


to find out which of the two models, blackbody or hy¬ 
drogen atmosphere, is more appropriate to describe the 
pulsar spectrum. If the future data will favour the at¬ 
mosphere model and will confirm that the surface tem¬ 
perature is as low as it follows from the current data, 
then it will result in important consequences for under¬ 
standing the physics of the neutrino emission in dense 
cores of NSs. 
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Figure 6. Top: ID and 2D marginal posterior distributions for A^h, T and D in the BB+PL model without the prior. Distance D is 
shown in units of The solid, dashed and dot-dashed lines in Nn-D frames show the empirical relations assuming R = 16 km, 

1 km and 20 km, respectively. See text for details. Vertical dashed lines in ID distributions show 5%, 50% and 95% quantiles. Bottom: 
The ID and 2D marginal posterior distributions for Vh, T and D in NSMAX+PL model without prior. Other options are the same as 
in the top panel. The temperature T is redshifted. 
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Figure 7. ID and 2D marginal posterior distributions for Nn, T, R and D in the BB+PL model (top) and the NSMAX+PL model 
(bottom), with account for the prior. Other options are the same as in Figure [6] 
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